On the stability of accelerating relativistic shock waves 



Giuseppe Palma and Mario Vietri 

Scuola Normale Superiore, Pisa, Italy 

ABSTRACT 

O 

^ . We consider the corrugation instability of the self-similar flow with an accel- 

erating shock in the highly relativistic regime. We derive the correct dispersion 
relation for the proper modes in the self-similar regime, and conclude that this 
solution is unstable. 
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1. Introduction 



It has long been known that the propagation of a shock through the outer layers of a 
star may, for sufficiently steep density distributions, lead to shock acceleration. Self-similar 
solutions (Gandel'man and Frank-Kamenetskii 1956, Sakurai 1960, Raizer 1964, Grover and 
Hardy 1966, Hayes 1968) notoriously provide excellent approximations to the generic solu- 
tions, for a reason most clearly discussed by Raizer (1964): behind the shock, a sonic point is 
formed, which prevents causal connection between the region immediately behind the shock 
with the more distant regions, which are sensitive to boundary conditions. Disconnection 
from initial conditions limits the number of parameters on which the solution may depend, 
while still allowing the existence of a self-similar solution. Because of the independence 
from initial conditions, it is expected that, shortly after their formation, all solutions with 
accelerating shocks will quickly approach the self-similar solution. 

Because of the possible relevance of these solutions to the evolution of Gamma Ray 
Bursts, there has been recently a strong surge of interest in their properties, especially in 
the relativistic regime, which, at least asymptotically, should be most appropriate to GRBs. 
A highly relativistic, self-similar solution has been presented by Perna and Vietri (2002) 
in planar geometry, for which accelerating, self-similar solutions appear for exponentially 
decreasing density stratifications outside the star. In spherical geometry, it is instead possible 
to find a richer spectrum of self-similar solutions, for basically any type of power-law density 
stratification (Best and Sari 2000, Pan and Sari 2005, Sari 2005). These solutions are often 
referred to as Type II solutions (Barenblatt 1996), meaning that the behaviour of the full 
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solution cannot be determined from dimensional analysis alone, like in the well-known Sedov 
solution, but is fixed instead by the conditions of regularity at some particularly difficult 
point, like the sonic point. 

The stability of relativistic hydrodynamic solutions with shocks is rarely considered in 
the literature. Gruzinov (2001) considered the linear stability of the spherically symmetric, 
relativistic explosions of Blandford and McKee (1976), but nothing is currently known about 
the stability of accelerating solutions, in the relativistic regime. In the Newtonian regime, 
the self-similar solutions are known to be unstable, from both a linear analysis (Chevalier 
1990), and a numerical, non-linear one (Luo and Chevalier 1994). The relativistic solutions 
differ from this analysis in two distinct respects: on the one hand, they include relativistic 
effects, on the other one, since the shock is supposed to be highly relativistic, they are 
obliged to assume the validity of the ultra-relativistic equation of state, p = e/3. Both of 
these circumstances lead to a greater susceptibility to corrugational instabilities. In fact, 
softening the equation of state leads to corrugational instability even for shocks propagating 
into a constant density environment: Ryu and Vishniac (1987) found a range of instabilities 
provided the adiabatic index is 7 < 1.2. Also, assuming an adiabatic index close to 7 = 5/3 
leads to both a narrower instability range in the parameter k, and in smaller growth rates 
than determined by Chevalier (1990) for 7 = 4/3 in the Newtonian, accelerating self-similar 
flow (Vietri, unpublished). 

On the other hand, relativistic effects tend to make pressure waves slow when compared 
to bulk motion, and thus to decouple nearby layer of fluid. Thus, if a purely kinematic 
difference exists in the motion of adjacent layers, it is more likely to grow because of the lack 
of restoring effects. 

The plan of this paper is as follows. In the next section, we shall briefly summarize the 
properties of the zero-th order solution of Perna and Vietri (2002). In the following section 
we shall perform a linear stability analysis, showing that accelerating, self-similar shocks are 
indeed unstable. In the last section, we shall discuss our results. 

2. The unperturbed solution 

Here we briefly recall the relevant properties of the zero-th order solution. The self- 
similar solution for an accelerating hyperrelativistic shock wave in a planar exponential 
atmosphere, in the highly relativistic regime, was derived by Perna and Vietri (2002). In 
their notation, the ambient density is given by 



p(x) = p e 



(1) 
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where po and ko are constants. If X l (s) denotes the shock position at proper time s, dimen- 
sional and covariance arguments impose 

dX l ak % , 

where = (0, /co, 0, 0) in the upstream frame and a is a dimensionless (negative) constant 
to be determined by imposing a smooth passage of the flow through a critical point. 

Integrating (2), it follows that 

— = log^ T - T7 j -- + const, (3) 



whose hyperrelativistic limit is 



-log— (4) 



-a ° Tj 



(V denotes the shock speed, T its Lorentz factor and the subscript i refers to the initial 
condition). Moreover, from (1), 

Ff ~^ . (5) 



Ti VPi 

In order to determine the value of a, the exact adiabatic fluid flow equations, 

(6) 

(nV)^ = 1 ; 

(y^ = wu^u u — pg^ v is the energy-momentum tensor, w and p being respectively local 
proper enthalpy density and pressure, is the fluid four-speed and n' is the baryon number 
density as measured in the comoving frame) (Landau and Lifshitz, 1989), are considered in 
their highly relativistic limit given by 

D e 3 

4 fe - + (Ve + ^j=0; (8) 

^ + V-(„v)=0; (9) 

here v and 7 are respectively the fluid speed and Lorentz factor, and e = 3p = 3u>/4 is local 
proper energy density, where use of the relativistic equation of state, p = e/3 has been made, 
as appropriate in the limit of highly relativistic motion. The operator D / Dt is the usual 
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convective derivative, and n = n lr y is the baryon number density as seen from the upstream 
frame, the comoving density being n' . 

The hyperrelativistic limit of Taub's jump conditions across a planar shock, with van- 
ishing speed parallel to the shock surface, is given by 

n 2 = 2r 2 m, 7 2 2 = ^r 2 , e 2 = 2rV, (10) 

where the subscripts 1 and 2 respectively refer to immediately pre- and post-shock quantities 
(Landau and Lifshitz, 1989). Upstream matter is obviously cold in this limit, since its sound 
speed obeys c s <C c; thus, using pi <C e± ~ p±, from (5) it follows that 

e 2 = 2^T 2+a = 2q T 2+a (11) 
i 



and 



n 2 = 2^r 2+ * = 2^ r 2+Q , (12) 



where n = po/m. 

Adopting as self-similarity variable 

f = k (x-X(t))T\t), (13) 

the equations (10), (11) and (12) can be cast into self-similar form by means of the definitions: 

7 2 (x, t) = g(0T 2 (t), e(x, t) = q R(OT 2+a (t), n(x, t) = z N(OT 2+a (t), (14) 

with 

9(0) = -, R(0) = 2, N(0) = 2. (15) 

Substituting (13) and (14) into (7), (8) and (9) it is possible to write equations for g, 
R and N in the form of a Cauchy problem with ordinary derivatives (the prime stands for 
the first derivative with respect to £): 

= 2g [-2a{4 + a ) + (2 + a)(a- 4Qg\ R 

a 2 + {a-^)g[-4a + {a-^)g] 1 J 

,_ g 2 [4{a-^)g-Ua-Za 2 } 



a 2 + (a-40g[-4a + (a-40g] 
= Ng K(18 + 5a) + 2(q - 4Qg [-2a(5 + 2a) + (2 + a) (a - 4Qg]} 

_ a 3 + ( a _^{ 5a 2 + ( a _ 4 ^[_ 5a + ( a _^]} \ ) 
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2.1. Analytical solution 

At the referee's request, we show here how to obtain a (nearly complete) analytical 
solution, which may perhaps prove useful. 

First of all, we determine a, by demanding the simultaneous vanishing of the numerators 
and of the denominator of eqs. 16 and 17. We find: 

a = -(2 + 4/V3) ~ -4.309401 . (19) 

Now we introduce a more convenient quantity, 

y = (a -40*7(0, (20) 

for which we know both the value at the shock front, Y , and at the critical point Y c (from 
the vanishing of the numerator of eq. 16, for instance): 

n = f, Y c = a(2-V3). (21) 

From eq. 17 we also derive the following differential equation for Y: 

Y (2a - 3a 2 )Y - 4a 2 
a-Ai a 2 + Y 2 — AaY ' 1 ' 

which can easily be integrated to yield the solution in implicit form: 

T r A f V & a 2 + Y 2 -4aY 
a-4e = aexp/ 1 , 1 = ~y (2a — 3a 2 )Y — Aa 2 ' (23) 

where suitable boundary conditions at the shock (£ = 0) have already been inserted. We 
also find: 

T f ( w ( /0 x f ( , 4 f y a 3a(-4 + 12a + 3a 2 ) ln(4a - 2y + 3ay) \ 
h = h(y)-h(a/2) , h(y) = - [^- 2 + I lny ) 

(24) 

If we now compute I\ for y = Y c , we can find an explicit expression for £ c , the location of 
the sonic point. After some labor, we find: 

1 1 g-6+7^/2 

& = ---_+ — » -0.462962... , (25) 

exactly the same value determined numerically in Perna and Vietri 2002. 

Now that we have found the dependence of Y (or of g, which amounts to the same) 
on £, we can use eqs. 16 and 18 to express R and as functions of Y, using the highly 
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convenient fact that Rl and N' do not depend upon £ except than through the combination 
Y. We obtain: 

dR _ (2 + a)Y-2a(A + a) 

dY (2a - 3a 2 )F - 4a 2 1 ' 

with the obvious solution (with the right boundary conditions at the shock already inserted): 

f Y (2 + a)y - 2a(4 + o 
/ a/2 <4/ (2a - 3a 2 )y - 4a 2 



f Y , (2 + a)y - 2aU + a) , x 

R = 2expl 2 , I 2 = 2 dy K — } \ K A J (27) 



from which we easily obtain 



R = 2 exp ^2(^ - 2^3)(y - a/2)^j 



(28) 



Proceeding analogously for N, we find: 

1 dN a 2 + Y(Y - 4a) a 2 (18 + 5a) + 2Y [-2a(5 + 2a) + (2 + a)Y] 
NdY ~ (2a-3a 2 )F-4a 2 -a 3 + Y [5a 2 + Y (-5a + F)] 

Including boundary conditions, this has the solution: 



(29) 



a/2 a(a-y)(-2y + a(A + 3y)) 



We also find: 



h = fs(Y) - h(a/2) , f 3 (Y) = _ (-2(2 + a)(-2 + 3a)F+ 

4a 2 (-4 - 4a + 57a 2 + 9a 3 ) arctan( 3a2 + 4y - 6a(1+y) ) 

yj — a 2 (2+3o) 2 

v/-a 2 (2 + 3ap + 
a(-8 + 22a + 9a 2 ) ln(-2F 2 + 3aF(2 + Y) - a 2 (4 + 3Y))) . (31) 

We remark that this is only a partially analytic solution in any case, because the rela- 
tionship between £ (the radial coordinate) and g or Y is implicit, making it easy to derive 
the value £ c of the sonic radius, but making it of little practical importance anyway. 

It is now possible to look at this solution in an instructive physical way. In fact, one 
may wonder how a solution may be accelerating, when a finite amount of energy is available: 
the answer is energy concentration. We can compute the parameter 

_ dlnE 

K = — — 32 
dlnM v ' 
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for, say, all the matter inside the sonic point. It is easy to check, thanks to the formulae 
given above, that E oc T 2+a , while M oc T a , both coefficients of proportionality being of 
course independent of time. Thus 

U = 1 + - « 0.535855 . (33) 
a 

From this we see that the specific energy (per baryon) increases as oc M-° Ae , and thus 
it increases without bounds when the total number of baryons between the shock and the 
sonic point decreases. At the same time, the total amount of energy in this layer of matter is 
dwindling to zero, since it is proportional to oc T 2+a = p -2 - 309 -: smaller and smaller amounts 
of total energy are used to propel even smaller amounts of matter. 



3. Perturbation analysis 

There are two ways to carry out the perturbation analysis. Here we present a simple 
approach, while in the Appendix we carry out the full job of perturbing both the equations 
of hydrodynamics and the Taub conditions at the shock, putting them together, and finding 
the only non-singular solution to the resulting set of equations. The reason for this double 
approach is that the first one is commendable in its simplicity, and is fully satisfactory for the 
aims of this paper, while the second one provides some crucial details which are necessary 
when comparing the early development of the numerical solutions to the fully non-linear 
problem with the linear, semi-analytic solution; the numerical solutions will be presented 
elsewhere, but we can anticipate that these details will come handy there. 

The argument runs as follows. It has been remarked before (Wang, Loeb and Waxman 
2003) that the perturbation problem is not self-similar, for the following reason. There is 
an intrinsic scale-length, l/fco, in the problem, which also determines the typical transverse 
length-scale on which perturbations in the shocked fluid can travel: this is l/(Tko), because 
the time-scale for a fluid element moving at speed 1 — r~ 2 /2 to cover the distance 1/fco 
(in the upstream frame) is l/(k T) in its frame; this implies that the maximum transverse 
length coverable by the fluid perturbation is l/(fc r) in both fluid and observer's frames. 
The ratio between this transverse length, and the transverse wavelength 1/k of the pertur- 
bations is adimensional, but does depend upon t through T: it thus breaks the well-known 
theorem according to which all adimensional quantities, in self-similar solutions, must be 
time-independent. 

Yet, we can consider two distinct regimes. The first one is that of short wavelenghts, 
i.e. k/ko 1. In this case the shock moves as if in a homogeneous medium, in which case it 
is well-known (Anile and Russo 1985) that no instability is present. Alternatively, we may 
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consider the opposite limit, k/(koT) <C 1, in which case a self-similar perturbation analysis 
is warranted: thus, we may consider the case k = as an approximation to those late times 
when, T being ^> 1, causal phenomena transverse to the shock's direction of motion cannot 
carry disturbances too far. In other words, we neglect the presence of causal phenomena 
between crests and valleys, and suppose them to evolve independently. 

Thus, each y = constant slice in the fluid evolves as an independent solution of the 
unperturbed equations, with slightly perturbed constants in the equation for the shock lo- 
cation: the equation we found above, eq. 2, can be integrated to give, to most significant 
orders in T: 

x = t -%(£)' + «' r = r,«p<-£). (34) 

We obtain distinct modes if we perturb the two constants, Tj and c±, between crests and 
valleys. When we perturb the quantity C\ in eq. 34, we are leaving the speed difference 
between crests and valleys unaffected, and thus SX, the shock displacement with respect to 
the unperturbed solution, must be independent of time. Alternatively, we may perturb the 
quantity Tj in eq. 34, in which case we find 5X oc T~ 2 . These are the two independent modes 
of the problem. 

One can likewise derive the expressions for the perturbed quantities behind the shock. 
Let us consider for instance & (identical calculations hold for Sn): we use 

e = 2T 2 p(X)R(0 . (35) 

For the first perturbation, we obviously have ST = 0, and thus <5£ = — k 5X r 2 , from which 
it follows that. 

& = 2T 2 (~k oP R5X + pR'(-k 5X T 2 )) « -2k 5X T 4 pR' . (36) 

This equation shows both the time dependence of 5e for this mode (oc T 4 p(X) oc T 4+a ), and 
its radial profile oc Rl (while 6n oc N'). We remark that, in this case, the fractional energy 
(density) perturbation grows as 

— = — oc T , (37) 

e n 

and it thus grows considerably. Analogously, 

7 2 = r 2 s(0; (38) 

V = rV(? = -koSXoT^g'. (39) 
For the other mode, 5X = 0, while ST/T = Sr /T . It follows that 
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(«) 

and, neglecting exponentially small terms, we find: 

5e = Uv 2 pR - pRa + 2T 2 pR' - ) ^ « r 2 p(4i? + 4£# - a#) ^ . (42) 

From this we again recover both the time dependence of the mode (oc T 2 p(X) oc T 2+a , which 
gives this time Se/e independent of time), and its space dependence, oc 4R + 4£i?' — aR' 
(identical equations, with the obvious substitutions, hold again for Sn). 

For the Lorentz factor, we find: 

V = r 2 ^(2, + 2^-^') (43) 



Su x ~ 5y 


oc 


r 1 


be 


oc 


p2+, 


5n 


oc 


p2+, 


5X 


oc 


p-2 



4. Discussion 

To summarize, we stress here that the hydrodynamical quantities grow as follows. Let 
us call u x the ^-component of four-velocity. For 5T ^ 0: 

(44) 
(45) 
(46) 
(47) 

while for the other, stronger mode, 8X ^ 0, we find: 

(48) 
(49) 
(50) 
(51) 

Obviously, the most surprising result is that, in this last mode, the concentrations of energy 
and baryon number be and Sn grow very fast, as oc r 4+a , or, stated otherwise, 

be bn ~ . . 

— oc — oc T 2 . 52 

e n 

Obviously, these concentrations will quickly become nonlinear, and their subsequent fate can 
only be ascertained through a numerical analysis. The same conclusion is reached when one 
considers the quick growth of bu x oc T 3 . 



5u x ~ by 


oc 


r 3 


be 


oc 


p4+, 


Sn 


oc 




5X 


oc 


r°. 
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Fig. 1. — Spatial dependence of <5y, when s = 3. In both this figure and in the following 
ones, we plot both the result of the analysis given in the text in their analytical form, and 
that given by the full analysis given in the Appendix: that the two superpose within the 
curve thickness bears witness to the complete equivalence of the two distinct approaches. 



The linear growth rates are so high that we can expect basically all small perturbations 
to become nonlinear during the acceleration of the ejecta of a Hypernova, when T changes 
from T i - 1< 1 to T f « 100. 

We should remark that the independence of these results from k is a peculiarity of 
the relativistic solution which does not exist in the Newtonian counterpart of the problem 
(Chevalier 1990). This independence is derived in an explicit way in the Appendix. 

Our conclusions differ from those of Wang, Loeb and Waxman (2003), who were unable 
to find self-similar perturbations to a similar problem in spherical geometry, reporting only a 
marginal growth for perturbations in the linear regime for an accelerating shock in spherical 
symmetry. Though not exactly identical, the growth rates should however coincide in the 
limit of large wavenumbers, where curvature effects become negligible. In support of our 
work we can make two different points. On the one hand, there is the the coincidence of the 
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Fig. 2. — Spatial dependence of Se, when s = 3. 

growth rates when these are computed in two, very different ways: the simple, physical one 
(presented above) and the more detailed, mathematical one to be presented in the Appendix. 

On the other hand, there is the similarity between the stability analysis in the New- 
tonian (Chevalier 1990) and relativistic regime (this paper): like us, Chevalier found two 
independent modes, for sufficiently small wavenumbers (the limit in which we can compare 
directly the two sets of computations), with very simple indices, s = 0, 1 (but please notice 
that they are defined in a different way from ours). Of these, the mode s = is of course 
only marginally unstable, while the other one, corresponding to a time growth oc for 
t — > 0, is identified by Chevalier (1990) as the physically relevant one. This is similar to 
our result of two independent modes, with distinct indices Se/e oc V s , s = 0, 2, the one with 
s = 2 resulting in the most severe (and thus physically relevant) instability. 

In a future paper, we will investigate the nonlinear development of the corrugational 
instability discussed here. 

Thanks are due to the referee, Dr. R. Sari, for helpful comments. 
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Fig. 3. — Spatial dependence of 5n, when s = 3. 



Appendix 

We present here, in a succinct form, the full perturbation analysis. 

A perturbation wrinkling the planar surface of the shock studied above causes, because of 
the refraction of the flow lines crossing an oblique shock, a transverse component (along y) in 
the shocked matter speed. We perform here a self-similar analysis of the small perturbations 
induced by the shock corrugation. In other words, we take, for the perturbed quantities: 

Su x (t, e, y) = u x (t, £, y) - u xo (t, = g x {i)V s {t)e ik ^ (53) 

5u y (t,£,y) = uy(t,£,y) = g y {i)Y s+ ^{t)e tk y (54) 

6e(t,£,y) = e(t,£,y) - e (t,0 = q Ri(iW s+a+P2 {t)e ik » (55) 

5n(t, t, y) = n{t, y) - n (t, = ^i(0r s+a+P3 (t)e^ . (56) 

s is the first-order analogous of the Oth-order parameter a and it will be determined by 
performing a new critical point analysis. The free parameters p±, p 2 and p 3 will be determined 
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Fig. 4. — Spatial dependence of 5y, when s — 1. 



shortly. In order to perturb the shock jump conditions, we will also need a self-similar form 
for the wrinkle (always in the upstream frame): 



5X(t,y) = X(t,y) - X (t) = -r' + w(*) e *w. 

K 



(57) 



Also ^4 will be determined shortly. To simplify the writing of the equation, we will define g 
which differs from Perna and Vietri's: 



fl'new = v^old- 



(58) 



We can also obtain linearized expressions for two useful quantities: in the hyperrela- 
tivistic limit, 



h(t,t,y) = $u x (t,£,y) 1 



1 



27o 2 (t,0 



8v(t,£,y) 



1 f5u x (t,£,y) 



7o(*,0 V 7 2 (U) 



x + 6uy(t,£,y)y 



(59) 
(60) 
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Fig. 5. — Spatial dependence of 5e, when s — 1. 



Now we linearize (7), (8) and (9): 



dSn d . . . . d . 

(8 7o e 57 + 4 7o 2 5e) + v ^) + 
47o 2 eo ( 9 ~t + 1&8v x + v 9 -t) + f + ^of + %6v x = 



47o 2 e 



9 9 
<9t dx 



d5v y 

4/3 



+ v - 



88v, 



dx 



85e de 



1/3 



3n, 



7/3 



n, 



4/3 



3n 



4/3 



+ 5v a 



<9 7o e 



0. 



(61) 

(62) 
(63) 

(64) 



Substitution of (53), (54), (55) and (56) into (61), (62), (63) and (64) yields (after some 
algebra and neglecting terms of manisfestly inferior order in T): 



p>s-i0 {-2aN ig ' + agN{ + g 3 [2(p 3 + s + a)N ± - (a - ^)N[}} 



-2a 



gg x N' + N (i±T^g 3 g y - 3g x g' + gg' x ) 



= 



(65) 
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Fig. 6. — Spatial dependence of 8n, when s — 1. 
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4R K [2(s - 1)£ 3 + ((a - 4£)<7 2 - 3a) <?'] - «? [-a + (a - 4£)<7 2 ] ^} 
+2agg x R' + T^g {2R l [(4 + p 2 + s + a)^ 3 - 2 (-a + (a - 4£)<? 2 ) </] 
- 5 [a + (a-4O^ 2 ] J R' 1 } = 



(66) 



A ar ^- V^i - r" 1 {2R [g y (2(p! + s - 1)(? 3 + 



-a 



-g (-a + (a- 4£)g 2 ) g' y ] + agg y R'} = 



2N {4g [a + 2(a - 40<? 2 ] fc^' + iV [2/2 (2(3s -a- l)g 3 
- (5a + (a - 40<? 2 ) <?') - Sg (-a + (a - 4£)<7 2 ) g' x ) 
-3g (a + 2(a-40g 2 )g x R'}} 
ST^Ng {2Ri [2g (a - (a - 4£)# 2 ) TV' + iV ((4 - 3p 2 - 3s + a)g 3 
-2ag' + 2(a - 40<?V)] + 3Ng [-a + (a - 4£)g 2 ] R[} 
+4P*- 1 {7# [a - (a - 4£)# 2 ] iVi/W + iV [i? (-2(3p 3 + 3s - a - 4)^ 3 iVi 
+ (4AW + 3^) (-a + (a - 4£)<7 2 )) + 3<? (-a + (a - 4£)<7 2 ) ^i?']} = 0. 



(67) 



(68) 
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We now determine Pi,P2,P3- We know from the non-relativistic problem (Landau and 
Lifshitz 1989) that the shock corrugation introduces three kinds of perturbations into the 
post-shock flow: entropy perturbations, vorticity perturbations, and sound waves. The first 
arise because, depending on the instantaneous state of motion of the corrugated surface, 
matter to be shocked may be faster or slower than the average flow; this in turn means 
that, after the shock, this matter may be hotter or cooler than average, and this implies 
entropy perturbations. Also, refraction of flow lines from oblique shocks produces vorticity 
perturbations. And lastly, pressure waves directed away from the shock are all but inevitable. 
The amplitude of these perturbations are coupled by the perturbed shock jump conditions, so 
that all physical quantities in the post-shock flow appear as the superposition of three kinds 
of perturbations. This simply implies that, in the above equations, we should choose the 
three coefficients Pi,P2,P3 hi such a way that no physical quantity is always negligible. Such 
situations are allowed by the above equations, but they describe perturbed flows where fewer 
independent perturbations are present. This of course may occur because the above equations 
still bear no information about the shock jump conditions, and thus may potentially refer 
also to strictly local perturbations. Inspection of the above equations shows that the only 
satisfactory solution where no physical quantity is completely neglected, as suggested by 
physical intuition, is: 



The fact that one solution for three parameters with four equations to be satisfied can be 
found bears witness to the soundness of the idea that the perturbed flow is also self-similar. 

Neglecting terms of lower order in T, it is possible to rewrite the set of differential 
equations as a 4D-Cauchy problem: 



x {7g 2 [-a 2 + (a- 4£)V] N X RN' -Ng[a + (a- 4£)g 2 ] 
[4(a - At)g 2 R{g x N> + A^') - 2aR(2g x N' + 5N l9 ') - 3ag(R 1 N' + N ± R') 
+g 3 (3(a - 4f + #1 (8(1 + a)R + 3(a - 4£)R'))] 
+N 2 [6(4 + a)(a- A^)g 6 Ri + Ua 2 g x Rg' - 48a(a - ^)g 2 g x Rg' 
+g 4 (-3a(7 + 3s + a)R 1 + 10(a - ^) 2 g x Rg') - 3a 2 g(AR l9 ' + g x Rl) 
+(a - A^g 5 (-6(a - 4£)Rrf + g x (2(6s + a - 8)R + 3(a - 4^^)) 
+ag 3 (18(a - ^)Rig' + g x ((20 - 24s + 2a)R + 9(a - 4£)R'))]} (70) 



P2 = P3 



= 1 . 



(69) 





2Rg v [2(s-3)g 3 



ag'+{a-iii)9 2 g l ]+agy9R'-ig' 2 ^ ) aR 1 



(71) 



2gR[- a +( a -^)g 2 \ 
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r>l _ 2 

^1 3Af 2 S 2 [a 2 -4o(a-4$)g 2 +(Q-45) 2 9 4 ] 

x j-7# 2 [a - (a - A()g 2 } 2 NxRN' - N 2 [3(3 - s + a) (a - 4£)g*Ri + 4a 2 g x Rg' 
-20a(a - Ai)g 2 g x Rg' + Qa 2 gR l9 ' + g 4 (6(1 + s)aR 1 + 4(a - ^^ifo') 
-2a# 3 ((4 - 6s + + 3(a - 4£)i? 1 # / ) 

+(a - ^)g b g x (2(a - 2)R + 3(a - 40^)1 + ^ h« + (« " 4£)<7 2 ] 
[4(a - 4£)g 2 R(g x N' + Aq«/) - 2aR(2g x N' + 5Aq«?') - ZagfaN' + N.R') 

+g 3 (3(a - 4^)i?iJV' + iVi (8(1 + a)i2 + 3(a - 4£)i2'))]} (72) 



/V' = 1 

JV 1 3Afg 2 [-a 3 +5a 2 (o-409 2 -5a(a-45) 2 3 4 +(a-4^) 3 g 6 ] 

x {Tag 2 [a 2 -{a- 4f) V] N^N' + aN 2 [-6(4 + a)(a - 4£)# 6 J Ri 
+4a 2 ^^' - 24a(« - 4£)g 2 g x Rg' + g 4 (3a(7 + 3s + 
+8(a - 4g) 2 g x Rg') + Zo?g{m lS f + 
+a# 3 (-18(a - 4£) J Ri£ / + & ((24s - 2a - 20)i? - 3(a - 4f)#))] 

+7V</ [2a 2 (7 + 3s + 7a)ff 3 JVii2 + 6(1 + s + a) (a - 4£) 2 g 7 N 1 R 
-2a(a - ^) 2 g 4 R{g x N> + A^') + 6a 2 (a - 4£)g*R(4g x N' + 3Aq<?') 
-2a s R(5g x N' + 8N ig ') - 3a 3 £(i2iiV + A^i?') 
+a(a - 4£)# 5 (3(a - 4£)#iA r/ + A^i (-8(2 + 3s + 2a)R + 3(a - 4f)#))]} ■ ( 73 ) 

These equations are linear in the perturbations, so their denominators are completely de- 
termined by the Oth-order solution. This means that it is possible to find roots of these 
terms without integrating (70), (71), (72) and (73). It is straightforward to show that 
the denominator of (71) never vanishes: since <7 2 (£) < 1/2 and a,£ < 0, it follows that 
-a + (a - 4£)y 2 > 0. 

We now remark that both the term common to the denominators of (70) and (72) 
and the denominator of (73) vanish at the Oth-order critical point £<> We have thus to 
impose that, by means of a judicious choice of a unique parameter s, the numerators of three 
equations vanish simultaneously with their denominators. That this can be done at all does 
appear a bit miraculous. 

We now turn to the perturbation of the shock jump (Taub) conditions, which will provide 
boundary conditions for the numerical integration of the above perturbed flow equations. 
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Exactly like in the Oth-order problem, Taub conditions across the shock provide bound- 
ary conditions 5^(0), g y (0), R\(0) and iVi(O). However, perturbations add some complica- 
tions. In primis, the shock does not have a unique speed v , but its points move along the 
x axis with y coordinate dependent velocities: 



v(y,t)=v (t) + SX(y,t)x. 



(74) 



Considering SX as a small perturbation, in the hyperrelativistic limit 1 , the perturbed Lorentz 
factor of the shock 2 is 



r up (y,i) = r(i) 



' -w :, / ( r 2 (0 - \ 



(75) 



Secondly, because of the wrinkle, the shock normal versor n does not coincide with x 
(such as the tangent versor t ^ y ). It is now possible to determine the boundary conditions 
arising from the shock jump conditions. We consider a frame locally (in the yc 3 coordinate) 
comoving with the shock, and denote quantities in this frame by means of primes; the usual 
conditions for the continuity of the fluxes of particle number, momentum and energy are: 
([/] means j\ - f 2 ) 

M = M = (76) 
= (77) 

= [wu' n u' t ] = (78) 

= [wiu' n ] = 0; (79) 

here n indicates baryonic density as measured in the frame locally comoving with the fluid; 
enthalpy density w and pressure p are always connected to energy density e in the usual 
hyperrelativistic way (Landau and Lifshitz 1989). Now we must connect these quantities 
with those for which equations (70), (71), (72) and (73) have been derived. 

Remembering the length dilation passing to the comoving frame, 

6X' yc (y,t) = r up (y c ,t)5X(y,t). (80) 

1 From now on, this approximation is always implicitly assumed. 

2 The subscript up reminds us that it is the same Lorentz factor with which one, in the shock frame, sees 
the upstream fluid incoming. 

3 Subscript c evocatively links to the comoving of the frame characterized by yc with the associated local 
shock surface element. 



P 
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If 5X< A;^ 1 , 



n(vc) 



9SX'y c (V,t) 



dy 



t(yc) 



dy 



V=VC 



(81) 



y=vc 



Beginning with the calculation of upstream quantities, spatial components of four-speed 
u\ = — r U pU are decomposed into normal and tangential parts: 



< = u\ ■ h = -r(l + 5XT 2 ) 
u'l = u\ ■ t = -r 



dSX' yc (y,t) 



dy 



(82) 
(83) 



y=vc 



Remembering that atmosphere is stratified, it follows that energy density is given by: 

ei = g r a (l-A; cDO, (84) 

while proper baryonic density is given by: 

ni = z a Y a {l - k 5X). (85) 



The Lorentz factor of the downstream fluid, as measured in the upstream frame, is: 



«2(f,S/c) = 7M) + 



dx 



&x(t,yc) + s>y(t,o,yc)- 



(86) 



x=0 



In a similar way, indicating with v^ o (t,x) the unperturbed downstream speed 1 — 2 ^ 2 | t x ^ , it 
is possible to write the spatial part of the four-speed as: 



u x 2 (t,y c )= 1 (t,0)v do (t,0) + 



d(l(t,x)v do (t,x)) 



dx 



5X(t,y c )+6u x (t,0,y c ) (87) 



x=0 



u y 2 (t,yc) = 5u y (t,0,y c ). 

Passing to the shock frame (as usual characterized by y c ), the four-speed transform as a 
four- vector: 



7 ^ M ' 2 = r up ( M °-M<) 

U*2 = r U p(«2 - \V\U° 2 ) 



Projecting u' 2 on n and i , it results 



u 2 y = u\. 



< = u' 2 x 



(89) 
(90) 
(91) 

(92) 



-20- 



d6X'(y,t) 
u« = u>* ' 



dy 



(93) 



y=vc 



Analogously, the energy density is 



e 2 = e (t, 0) + 



de (t, x) 



dx 



SX(t,y c ) + 5e(t,0,y c ), 



(94) 



x=0 



while the proper baryonic density (it is necessary to divide the density in the upstream frame 
by the downstream Lorentz factor u9,(t,yc)) is given by 



n o (t,0)+ 



dx 



n 2 = 



x=0 



SX(t,yc) + 6n(t,0,y c ) 



(95) 



Now it is possible to substitute the self-similar form for perturbations into (76), (77), 
(78) and (79). Linearizing Taub's conditions and neglecting terms of manifestly lower order 
in T, we find (after much algebra) 

4^2 p 4 -s + 2( v / 3 - 2)a] e - (l6g x (0) - V^N^ofj aT P4+3 = (96) 
4(2^3 - 3)e + 3 (±V2g x (0) - i?i(0)) P> 4+3 = (97) 



ike + V2k g y (0)r P4+s = 



2(s-p 4 ) + (11-6^) 



a 



e + 



(20v^^(0)-3i2i(0)) 



ar P4+3 = ^ 



(98) 
(99) 



The requirement that (96), (97), (98) and (99) form a non-singular system of equations for 
the boundary conditions imposes p 4 = — 3. Solving the system we obtain: 



9x(0) 



11 + 6^- (3 + 2^)5 
2^(7 + 4^) 

. , ike 
9y(0) - 



Ri(0) 

iVi(O) = 



V2~k 

2 [-27 - Uy/3 + (9 + 6-s/3)s] 
" 3(7 + 4^3) 

2 [-27 - 16>/3 + (21 + 12>/3)s] 
(3 + 2^(7 + 4^3) 



e. 



(100) 

(101) 
(102) 
(103) 
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All the perturbed variables scale with e: we are free to to set it = 1 in the following. 

The boundary conditions for g x , Ri and Aq do not depend upon k, nor do equations 
(70), (72), (73). Thus, the only quantity that does depend upon k is g y , but since, as 
remarked above, its denominator never vanishes, it follows that s cannot be fixed by g y nor, 
by implication, by k. Thus s is independent of k. 

We can thus restrict our discussion to the 3D Cauchy problem for g X) R x and Aq and 
look for a value of s which leads to 

lim N gx (0 = lim N Rl (£) = lim N Nl (£) = 0. (104) 

*-< *~< 

Here N gx means the numerator of the differential equation for g x , and likewise for N Rl and 
ATjvi- Boundary conditions are given in £ = 0, and it is from here that integration process 
can start leftwards, until it reaches the critical point £c < from the right. 

Using a binary search based on the direct applications of the theorem of zeroes, we 
studied the numerators near £ c varying s; we found these numerators could never vanish for 
complex s. Thus we investigated these solutions for real values of s. We find two solutions: 

si = 1 ; s 2 = 3 . (105) 

Integration of the system of equations for £ < £ c shows indeed that no divergence is present 
in our system of equations: see figures 7 and 8. 

It is also possible to insert the solutions derived in Section 3 into these equations, in 
order to check that the two sets of computations are mutually compatible. This has been 
done by means of Mathematica, since the computations are very heavy, and the expected 
mutual agreement has indeed been found. 

The dependence of all quantities upon the adimensional radius £ is shown in the figures, 
first for s = 3, and then for s — 1 (g y figures are obviously reported here as a result of this 
only section). It can easily be checked that we have indeed found two self-similar, distinct 
solutions, passing without divergence through the sonic point. 

Since s is real, so are g x , Ri and Aq. However, g y , the only function depending on k/ko, 
is purely imaginary: in fact, the boundary condition is purely imaginary and its differential 
equation contains only real terms multiplied by g y and a term containing R\ multiplied by 
i. Furthermore, ik/k enters g y (£) simply as a multiplicative factor. In a more intuitive way, 
let us consider several wrinkles of the same amplitude e/k , but with different wave-lengths. 
Clearly, the tangential component of the speed at the shock (which, after all, is the quantity 
that determines g y ) is proportional to sinarctane/c/A;o; as a wrinkle is a small perturbation, 
this factor reduces to ek/k , thus justifying previous mathematical result. 
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Fig. 7. — The run of numerators and denominators of g' x , R[ ed N[ around £c for s = 1. 
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